Competing roughening mechanisms in strained heteroepitaxy: a fast kinetic Monte 

Carlo study 
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We study the morphological evolution of strained heteroepitaxial films using kinetic Monte Carlo 
simulations in two dimensions. A novel Green's function approach, analogous to boundary integral 
methods, is used to calculate elastic energies efficiently. We observe island formation at low lattice 
misfit and high temperature that is consistent with the Asaro-Tiller-Grinfeld instability theory. At 
high misfit and low temperature, islands or pits form according to the nucleation theory of Tersoff 
and LeGoues. 
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Coherent three-dimensional (3D) islands in strained 
heteroepitaxial films are of great interest because they 
can self-assemble as quantum dots for possible advanced 
optoelectronic applications [Qj2|. They are observed in 
a variety of film-substrate combinations including Ge/Si, 
InAs/GaAs, InAs/InP, etc. In these systems, island for- 
mation follows the Stranski-Krastanov mode. Initially, 
two-dimensional (2D) layer-by-layer growth leads to a flat 
wetting layer under stress. Beyond a threshold film thick- 
ness, 3D islands emerge on top of the wetting layer par- 
tially relieving the stress. The precise island formation 
mechanism is currently under intensive debate. Accord- 
ing to the nucleation theory of Tersoff and LeGoues, the 
growth of stable islands requires overcoming an energy 
barrier associated with a critical island size How- 
ever, experiments reveal gradual development of ripples 
H or pre-pyramids || at the initial stage of island for- 
mation. These observations are more consistent with 
the Asaro-Tiller-Grinfeld (ATG) linear instability theory, 
which predicts that morphological perturbations at suffi- 
ciently long wavelengths grow spontaneously and steadily 
HQ. It was originally proposed for smooth surfaces, but 
extensions to faceted ones based on non-equilibrium de- 
position conditions or finite vicinality of substrates |[J 
have been suggested. 

To better understand the roughening mechanism of 
strained layers, we have performed kinetic Monte Carlo 
simulations using an atomistic model fl0|-|l3||. This ap- 
proach is computationally very intensive but can reliably 
account for both lattice discreteness and non-equilibrium 
conditions. Previous simulations successfully demon- 
strated island formation in strained layers but only via 
the nucleation mechanism |l0| - [T2]| . In this work we intro- 
duce significantly more efficient algorithms. Thus we can 
explore a much wider range of conditions and observe 
a rich variety of morphologies in better general agree- 
ment with experiments. In particular, by lowering the 
lattice misfit and raising the temperature, the roughen- 
ing mechanism crosses over from nucleation to instability 



controlled. 

We adopt the 2D ball and spring model of heteroepi- 
taxy defined on a square lattice first studied by Orr et 
al |n| and subsequently by Barabasi [[ll] and Khor and 
Das Sarma Simulations in 3D limited to submono- 
layer coverage were performed by Meixner et al fl3[| . Our 
model parameters are appropriate to the widely stud- 
ied Sii-zCe^/Si system. We assume a substrate lat- 
tice constant a s — 2. 71 5 A so that gives the correct 
atomic volume in crystalline silicon. The lattice constant 
cif of the film material is related to the lattice misfit 
e — (af — a s )/cif which has a compositional dependence 
e = 0.04x. Nearest and next nearest neighboring atoms 
are directly connected by elastic springs with force con- 
stants fcjv = 13.85eV/a^ and kpjN = kpj/2 respectively. 
This choice gives the correct modulus c\\ of silicon and 
a shear modulus constant along tangential and diago- 
nal directions, despite a slight anisotropy in the Young's 
modulus. The elastic couplings of adatoms with the rest 
of the system are weak and are completely neglected for 
better computational efficiency. Solid-on-solid conditions 
and atomic steps limited to at most two atoms high are 
assumed. Every topmost atom in the film can hop to a 
random topmost site s columns away where s = ±1, ±2, 
... or ±s mQX with equal probability. Previous simulations 
allowed only nearest neighbor hopping (i.e. s m ax = 1) 
pp| ILsj . To speed up the simulations, we put s max = 8 
or 20 respectively for x > 0.6 or x < 0.6. These hopping 
ranges are much shorter than the dimensions of the rele- 
vant structures (islands or pits) on the films and we have 
checked that decreasing s max does not alter our results. 
The hopping rate T m of a topmost atom m follows an 
Arrhcnius form: 



F m = R exp 



n m j - AE m - E 
k R T 



(1) 



Here, n m is the number of nearest and next nearest neigh- 
bors of atom m. We choose a bond strength 7 = OAeV 
which will be explained later. The energy AE m is the 
difference in the strain energy E s of the whole lattice at 
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mechanical equilibrium when the site is occupied ver- 
sus unoccupied. Finally, we put Eq = 0.53eV and 
Ro = 2D /(a s a s ) 2 with D = 5.2 x lO 13 !^ 1 and 
o\ = \{svnax + l)(2s maa; + 1). This gives the appro- 
priate adatom diffusion coefficient for silicon (100) |l4) . 
Our model follows detailed balance. 

The simulations involve intensive computations result- 
ing solely from the long-range nature of elastic interac- 
tions. Practically all the CPU time is spent on the re- 
peated calculations of E s which is needed to find AE m 
and hence T m in Eq. (Q). The elastic problem is for- 
mulated as follows. First, a flat film is homogeneously 
strained [Q. This provides a convenient reference posi- 
tion with displacement Hi = for every atom i. In gen- 
eral, the elastic force on atom i by a directly connected 
neighbor j is = — Ky (Ui — Uj) + bij after linearization 

ki-jhi-jn*: is the 



arises from the 



where the 2x2 symmetric matrix Ky 

modulus tensor and bij = — Zy)Kjjnij 
homogeneous stress in flat films. The spring constant fey 
equals either fcjv or n for tangential or diagonal con- 
nection respectively. The unit column vector hij points 
from the unstrained lattice position of atom j towards 
that of atom i and t denotes transpose. Furthermore, l®j 
and Uj are respectively the natural and homogeneously 
strained spring lengths which follow easily from a s and 
e. Mechanical equilibrium requires yU fij = for each 
atom i. This leads to a large set of equations coupling 
the Hi of all of the atoms. The solution then gives the 
elastic energy stored in every spring and hence E s . 

We now introduce a Green's function approach for cal- 
culating E s efficiently requiring the explicit considera- 
tion of only the surface atoms. It is a lattice analogue of 
boundary integral methods and is superior to boundary 
element techniques for our intrinsically discrete problem. 
We first derive the exact formalism. Figure 1 shows an 
example of a small lattice of atoms (solid circles). As a 
mathematical construct, we extend the lattice by adding 
ghost atoms (open circles) with similar elastic proper- 
ties. Unphysical couplings are hence introduced but can 
be exactly cancelled by applying external forces and 

f% to every real surface atom j and ghost surface atom 
j respectively with 



if = E( K 



33' u 3 



b J3 >) 



33 u 3 



(2) 
(3) 



The summation in Eq. (g) is over each ghost atom j' 
connected directly to the real atom j and it is analogous 
in Eq. (||). It is easy to see that the real atoms are then 
exactly decoupled from the ghost atoms which are held 
precisely at their homogeneously strained positions. 
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FIG. 1. Ghost atoms (open circles) are added on top 
of the real atoms (solid circles) forming an extended lattice. 
Unphysical couplings are exactly cancelled by external forces 
(arrows) applied to the real and ghost surface atoms. 

A lattice Green's function is then applied to express 
the displacement, Ui, of every real surface atom i under 
the influence of the external forces: 



5> 

j' 



(4) 



Note that the Green's function G is defined for the ex- 
tended lattice and is independent of the film morphology. 
It can thus be computed numerically once prior to the 
start of the simulation. It should not be confused with 
the half-plane Green's function which provides simpler 
but only approximate results |3 13 1. Combining Eqs. 
0) , we arrive at the reduced set of equations 



t, — / )[(Gjj - G,UK,,m/, - Gijbjf] 



(5) 



33' 



coupling only the real surface atoms where the sum is 
over all pairs of directly connected real and ghost surface 
atoms j and j' respectively. The solution of Eq. (||) gives 
u l . 

The elastic energy E s can then be calculated directly 
from 



(6) 



which is derived from a simple consideration of virtual 
work. Here, the sum is defined similarly to that in Eq. 
(||) and E^ equals the unrelaxed value of E s when u j = 
which can be straightforwardly computed. The method 
summarized in Eqs. (||) and (^) is exact and practical for 
simulations at moderate scales. 

We can go further and use a coarse-grained version 
of our Green's function approach for a further boost on 
the computational efficiency. Finding the strain energy 
AE m of the atom m needed in Eq. (|l|) requires cal- 
culating the strain energy E s of the whole lattice twice 
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with and without the atom m. Certain fine details of the 
surface far away are obviously unimportant and can be 
neglected. Specifically, surface atoms are grouped into 
sets with the ith of them denoted by fi/. We neglect 
fluctuations within a set by assuming identical displace- 
ment Hi = uj for each member i 6 fi/. Equation (||) is 
then approximated as 



Ui 



E 



E ( 



'1 3 



uj 



E« 



l3°33' 



(7) 



where Gij — Gij with the lattice point i at the cen- 
troid of the set Qi. Every atom within 3 columns from 
the atom m is not coarsened and constitutes their own 
single-membered set. Farther away at r columns from the 
atom m, sets contain atoms in neighborhoods of 2r/3 + 1 
columns wide, a form motivated by simple error analy- 
sis. We have checked numerically that a smaller degree of 
coarsening leads to no noticeable difference in our results. 

Surface diffusion can be simulated using the hopping 
rates in Eq. (|l|) as AE m is now readily computable. We 
adopt an acceptance-rejection algorithm aided by tabu- 
lated values of AE m for 5 8 sample surface configurations. 
Details will be explained elsewhere. In our main simula- 
tions, the lattice is of 1024 atoms wide following periodic 
boundary conditions. The substrate is 1024 monolayers 
(ML) thick while the extended lattice on which we com- 
pute the Green's function includes also a film of 80ML. 
Fixed boundary conditions are applied to the top and 
bottom layers of the extended lattice. 

We first simulate deposition of pure Ge film (i.e. x = 1) 
with misfit e = 4% at temperature T = 600K. At very 
high deposition rate R = 80MLs _1 [Fig. 2(a)], we ob- 
serve layer-by-layer growth. At slower deposition rate 
R = 8MLs _1 [Fig. 2(b)], the film is initially flat but pits 
then develop. A detailed examination of the morpholog- 
ical evolution indicates that the pits appear rather inde- 
pendently and suddenly. Once created, they are immedi- 
ately bounded by side-walls at an energetically favorable 
45° inclination. These features strongly support the nu- 
cleation mechanism for their formation, noting that pits 
are energetically more favorable than islands 0. This 
pit nucleation process is similar to that in Fig. 1(d) of 
Ref. go|. At R = 0.8MLS" 1 [Fig. 2(c)], islands with 45° 
side-walls nucleate at very early stage before the film is 
sufficiently thick for pit formation. The result is anal- 
ogous to those in Fig 1(a) of Ref. [|l0) and also Refs. 
p| , ^2[ . Further decreasing R towards realistic values of 
order O.OlMLs -1 leads to similar but more widely sepa- 
rated islands. 
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FIG. 2. Simulations of deposition of Ge films at T = 600^ 
with substrate width and thickness both equal 1024a s ~ 
2780A. The axes are in unit of a a 

Figure 3 shows results for deposition at misfit e = 2% 
with x = 0.5 at T = 1000/G Depending on R, we ob- 
serve analogous layer-by-layer growth [Fig. 3(a)], layer- 
by-layer growth followed by roughening [Fig. 3(b)], and 
island growth [Fig. 3(c)]. However, the islands in Figs. 
3(b)-(c) emerge gradually from ripple- like perturbations 
with local surface inclinations increasing steadily and 
relatively synchronously in agreement with experiments 
IimH) and ATG instability theory §,§. This regime 
has not been reported in previous atomistic simulations 
mainly due to inaccuracies in accounting for the long- 
range parts of the elastic interactions or the rather thin 
substrates used [ [To| pi) . Instead, it was observed in con- 
tinuum simulations ]17| , which however cannot realize 
the nucleation mechanism. 
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FIG. 3. Simulations of deposition of Sio.sGeo.s films at 
T = lOOOiG 

The importance of the lattice misfit in deciding the 
roughening mechanism is particularly easy to under- 
stand. The nucleation of islands or pits occurs at a rate 
Rnud ~ exp(— ce~ 4 ) with c being a constant Q and be- 
comes very slow at small e |l(| . The ATG instability with 
a roughening rate Rinst ~ e 8 @] then dominates. For 
Figs. 2(b) and 3(b), we have chosen deposition rates close 
to the relevant roughening rates. We then observe kinet- 
ically limited wetting layers prior to roughening which 
is characteristic of Stranski-Krastanov growth. However, 
the threshold thickness depends strongly on R contrary 
to experimental findings A more realistic model 

in the future should include other mechanisms such as 
film-substrate interactions |ll| or nonlinear elasticity |]] 
which have been argued to give a more stable wetting 
layer. 

We have also simulated annealing of initially flat films 
of 30ML at T = lOOOif . At this high temperature, rough- 
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ening is mainly due to the ATG instability, and we ob- 
serve the development of ripples followed by islands. Fig- 
ure 4 shows the surface profiles after the islands are fully 
developed. The cusps on the surfaces |ll|] are limited by 
either the substrate or the local step height limit imposed 
in our model. We have measured the island size I from 
power spectra of 11 realizations of similar surfaces. Wc 
obtain I ~ x -1 ' 8 in reasonable agreement with x~ 2 from 
the instability theory J6| but slightly different from x~ x 
from experiments Jl5| , ^6| . The discrepancy with experi- 
ments may be improved if the compositional dependence 
of film properties such as bond energies and diffusion co- 
efficients are properly considered. The island sizes in gen- 
eral lie within the experimental range due to our choice 
of the bond strength 7 = OAeV, which in fact is a rea- 
sonable value. 



ing description of island formation in Sii-^Ge^ films at 
low, and probably also at high lattice misfit How- 
ever, the nucleation mechanism applied to high misfit 
regimes in certain experimental situations has not been 
ruled out pfl| . Thus the competition of mechanisms can 
be important for interpreting experimental results. In 
our simulations, for deposition rates close to the rele- 
vant roughening rate, kinetically limited wetting layers 
develop prior to roughening. At lower but more realistic 
deposition rates, islands form at an early stage and are 
more widely separated. This may be related to the great 
variation in how closely the islands are packed under var- 
ious conditions in experiments 0j5|. 

We thank B.G. Orr for helpful comments. CHL is sup- 
ported by PolyU grant no. 5309/01P. 
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FIG. 4. Simulations of annealing of initially flat Sii-mGe^ 
films of 30ML at T = WOOK for a period of time t. 

We have already discussed possible extensions of the 
model to include film-substrate chemical interactions, 
nonlinear elasticity, and composition-dependent material 
properties. To further enhance the morphological resem- 
blance with experiments, one can explore more sophisti- 
cated forms of bond energies favoring 2D analogs of (105) 
and (113) facets {§. 

We should note that the validity of ATG instability 
theory in 3D is not clear due to a singular form of the 
equilibrium surface energy in the presence of facets 
Critical tests of the theory by 3D simulations are impor- 
tant, and may be feasible using our method. For exam- 
ple, our 2D simulation leading to Fig. 2(b) involves about 
6 x 10 7 hops and ran for 18 hours on a pentium 2GHz 
computer. Thus, extensions to 3D should be practical. 

In conclusion, using accelerated algorithms which 
properly and efficiently account for long-range elastic in- 
teractions, we have simulated deposition and annealing of 
strained heteroepitaxial layers in 2D. At low misfit and 
high temperature, we observe ripples and subsequently 
gradual island formation in consistent with the ATG in- 
stability theory. At high misfit and low temperature, 
islands or pits are generated via the nucleation pathway. 
These suggest a non-trivial competition between rough- 
ening mechanisms, although reliable quantitative deter- 
mination of the crossover conditions is beyond the scope 
of our model. The ATG instability is the most promis- 
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